rm(list=ls())

data <- read.csv("C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/Replication/masterdata_cces.csv")
attach(data)


## Table 1: Opinion across relevant demographic categories

agecat1 <- weighted.mean(incpolice[agecat==1],w=comwt[agecat==1],na.rm=T)
agecat2 <- weighted.mean(incpolice[agecat==2],w=comwt[agecat==2],na.rm=T)
agecat3 <- weighted.mean(incpolice[agecat==3],w=comwt[agecat==3],na.rm=T)
agecat4 <- weighted.mean(incpolice[agecat==4],w=comwt[agecat==4],na.rm=T)
blackmean <- weighted.mean(incpolice[black==1],w=comwt[black==1],na.rm=T)
whitemean <- weighted.mean(incpolice[white==1],w=comwt[white==1],na.rm=T)
hispmean <- weighted.mean(incpolice[hisp==1],w=comwt[hisp==1],na.rm=T)
asimean <- weighted.mean(incpolice[asi==1],w=comwt[asi==1],na.rm=T)
femalemean <- weighted.mean(incpolice[male==0],w=comwt[male==0],na.rm=T)
malemean <- weighted.mean(incpolice[male==1],w=comwt[male==1],na.rm=T)
nohs.mean <- weighted.mean(incpolice[nohs==1],w=comwt[nohs==1],na.rm=T)
hs.mean <- weighted.mean(incpolice[hs==1],w=comwt[hs==1],na.rm=T)
somecoll.mean <- weighted.mean(incpolice[somecoll==1],w=comwt[somecoll==1],na.rm=T)
coll.mean <- weighted.mean(incpolice[coll==1],w=comwt[coll==1],na.rm=T)
postgrad.mean <- weighted.mean(incpolice[postgrad==1],w=comwt[postgrad==1],na.rm=T)
dem.mean <- weighted.mean(incpolice[dem==1],w=comwt[dem==1],na.rm=T)
rep.mean <- weighted.mean(incpolice[repub==1],w=comwt[repub==1],na.rm=T)
indep.mean <- weighted.mean(incpolice[dem==0&repub==0],w=comwt[dem==0&repub==0],na.rm=T)
incunder30.mean <- weighted.mean(incpolice[incunder30==1],w=comwt[incunder30==1],na.rm=T)
inc30to60.mean <- weighted.mean(incpolice[inc30to60==1],w=comwt[inc30to60==1],na.rm=T)
inc60to100.mean <- weighted.mean(incpolice[inc60to100==1],w=comwt[inc60to100==1],na.rm=T)
incover100.mean <- weighted.mean(incpolice[incover100==1],w=comwt[incover100==1],na.rm=T)
urban.mean <- weighted.mean(incpolice[urban==1],w=comwt[urban==1],na.rm=T)
smallmet.mean <- weighted.mean(incpolice[smallmet==1],w=comwt[smallmet==1],na.rm=T)
rural.mean <- weighted.mean(incpolice[rural==1],w=comwt[rural==1],na.rm=T)

require(xtable)
xtable(as.data.frame(c(agecat1,agecat2,agecat3,agecat4,blackmean,whitemean,hispmean,asimean,
femalemean,malemean,nohs.mean,hs.mean,somecoll.mean,coll.mean,postgrad.mean,
incunder30.mean,inc30to60.mean,inc60to100.mean,incover100.mean)))

## FIGURE 1: Opinion by race and age

## SUMMARY STATS

black.agecat1 <- weighted.mean(incpolice[black==1&agecat==1],w=comwt[black==1&agecat==1],na.rm=T)
black.agecat1
black.agecat2 <- weighted.mean(incpolice[black==1&agecat==2],w=comwt[black==1&agecat==2],na.rm=T)
black.agecat2
black.agecat3 <- weighted.mean(incpolice[black==1&agecat==3],w=comwt[black==1&agecat==3],na.rm=T)
black.agecat3
black.agecat4 <- weighted.mean(incpolice[black==1&agecat==4],w=comwt[black==1&agecat==4],na.rm=T)
black.agecat4

black.agecat1.var <- (black.agecat1*(1-black.agecat1))/length(incpolice[black==1&agecat==1])
black.agecat1.var
black.agecat2.var <- (black.agecat2*(1-black.agecat2))/length(incpolice[black==1&agecat==2])
black.agecat2.var
black.agecat3.var <- (black.agecat3*(1-black.agecat3))/length(incpolice[black==1&agecat==3])
black.agecat3.var
black.agecat4.var <- (black.agecat4*(1-black.agecat4))/length(incpolice[black==1&agecat==4])
black.agecat4.var

white.agecat1 <- weighted.mean(incpolice[white==1&agecat==1],w=comwt[white==1&agecat==1],na.rm=T)
white.agecat1
white.agecat2 <- weighted.mean(incpolice[white==1&agecat==2],w=comwt[white==1&agecat==2],na.rm=T)
white.agecat2
white.agecat3 <- weighted.mean(incpolice[white==1&agecat==3],w=comwt[white==1&agecat==3],na.rm=T)
white.agecat3
white.agecat4 <- weighted.mean(incpolice[white==1&agecat==4],w=comwt[white==1&agecat==4],na.rm=T)
white.agecat4

white.agecat1.var <- (white.agecat1*(1-white.agecat1))/length(incpolice[white==1&agecat==1])
white.agecat1.var
white.agecat2.var <- (white.agecat2*(1-white.agecat2))/length(incpolice[white==1&agecat==2])
white.agecat2.var
white.agecat3.var <- (white.agecat3*(1-white.agecat3))/length(incpolice[white==1&agecat==3])
white.agecat3.var
white.agecat4.var <- (white.agecat4*(1-white.agecat4))/length(incpolice[white==1&agecat==4])
white.agecat4.var

hisp.agecat1 <- weighted.mean(incpolice[hisp==1&agecat==1],w=comwt[hisp==1&agecat==1],na.rm=T)
hisp.agecat1
hisp.agecat2 <- weighted.mean(incpolice[hisp==1&agecat==2],w=comwt[hisp==1&agecat==2],na.rm=T)
hisp.agecat2
hisp.agecat3 <- weighted.mean(incpolice[hisp==1&agecat==3],w=comwt[hisp==1&agecat==3],na.rm=T)
hisp.agecat3
hisp.agecat4 <- weighted.mean(incpolice[hisp==1&agecat==4],w=comwt[hisp==1&agecat==4],na.rm=T)
hisp.agecat4

hisp.agecat1.var <- (hisp.agecat1*(1-hisp.agecat1))/length(incpolice[hisp==1&agecat==1])
hisp.agecat1.var
hisp.agecat2.var <- (hisp.agecat2*(1-hisp.agecat2))/length(incpolice[hisp==1&agecat==2])
hisp.agecat2.var
hisp.agecat3.var <- (hisp.agecat3*(1-hisp.agecat3))/length(incpolice[hisp==1&agecat==3])
hisp.agecat3.var
hisp.agecat4.var <- (hisp.agecat4*(1-hisp.agecat4))/length(incpolice[hisp==1&agecat==4])
hisp.agecat4.var


asi.agecat1 <- weighted.mean(incpolice[asi==1&agecat==1],w=comwt[asi==1&agecat==1],na.rm=T)
asi.agecat1
asi.agecat2 <- weighted.mean(incpolice[asi==1&agecat==2],w=comwt[asi==1&agecat==2],na.rm=T)
asi.agecat2
asi.agecat3 <- weighted.mean(incpolice[asi==1&agecat==3],w=comwt[asi==1&agecat==3],na.rm=T)
asi.agecat3
asi.agecat4 <- weighted.mean(incpolice[asi==1&agecat==4],w=comwt[asi==1&agecat==4],na.rm=T)
asi.agecat4

asi.agecat1.var <- (asi.agecat1*(1-asi.agecat1))/length(incpolice[asi==1&agecat==1])
asi.agecat1.var
asi.agecat2.var <- (asi.agecat2*(1-asi.agecat2))/length(incpolice[asi==1&agecat==2])
asi.agecat2.var
asi.agecat3.var <- (asi.agecat3*(1-asi.agecat3))/length(incpolice[asi==1&agecat==3])
asi.agecat3.var
asi.agecat4.var <- (asi.agecat4*(1-asi.agecat4))/length(incpolice[asi==1&agecat==4])
asi.agecat4.var

black.opinion <- round(c(black.agecat1,black.agecat2,black.agecat3,black.agecat4),2)
white.opinion <- round(c(white.agecat1,white.agecat2,white.agecat3,white.agecat4),2)
hisp.opinion <- round(c(hisp.agecat1,hisp.agecat2,hisp.agecat3,hisp.agecat4),2)
asi.opinion <- round(c(asi.agecat1,asi.agecat2,asi.agecat3,asi.agecat4),2)

opinion <- c(white.opinion,black.opinion,hisp.opinion,asi.opinion)

Races <-c(rep("White", 4), rep("Black", 4), rep("Hispanic",4), rep("Asian",4))
agecats <- c("<30","30-44","45-64","65+","<30","30-44","45-64","65+","<30","30-44","45-64","65+","<30","30-44","45-64","65+")

dat <- data.frame(agecats,Races,opinion)
require(ggplot2)
p <- ggplot(dat,aes(agecats,opinion))
p +guides(fill=guide_legend(title=""))+xlab("") + ylab("")+ylim(0,.8)+
geom_bar(stat = "identity", aes(fill = Races),position="dodge")+
theme(
 panel.background = element_rect(fill = "gray85"), 
                 #  plot.background = element_rect(fill = "white", color = "white"), 
                  # panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                panel.grid.minor = element_blank(), 
                axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_text(size = 15), 
                 axis.text.y = element_text(size = 15)) +
          theme(
 panel.background = element_rect(fill = "white"), 
                 plot.background = element_rect(fill = "white", color = "white"), 
                panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                 panel.grid.minor = element_blank(), 
                  axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_text(size = 15), 
                  axis.text.y = element_text(size = 15),
                  axis.line.y = element_line(color = "gray40"),
          #   legend.position = c(0.9, 0.2),
                  plot.title = element_text(hjust = 0.2))

## FIGURE 2: Age-adjusted opinion among Asians, Blacks, and Hispanics

require(plyr)
shagecat1.white <- count(white[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat1.white
shagecat2.white <- count(white[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat2.white
shagecat3.white <- count(white[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat3.white
shagecat4.white <- count(white[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat4.white


shagecat1.black <- count(black[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat1.black
shagecat2.black<- count(black[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat2.black
shagecat3.black<- count(black[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat3.black
shagecat4.black<- count(black[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat4.black

shagecat1.hisp <- count(hisp[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat1.hisp
shagecat2.hisp<- count(hisp[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat2.hisp
shagecat3.hisp<- count(hisp[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat3.hisp
shagecat4.hisp<- count(hisp[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat4.hisp

shagecat1.asi <- count(asi[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat1.asi
shagecat2.asi<- count(asi[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat2.asi
shagecat3.asi<- count(asi[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat3.asi
shagecat4.asi<- count(asi[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat4.asi

shagecats.white <- c(shagecat1.white,shagecat2.white,shagecat3.white,shagecat4.white)
shagecats.black <- c(shagecat1.black,shagecat2.black,shagecat3.black,shagecat4.black)
shagecats.hisp<- c(shagecat1.hisp,shagecat2.hisp,shagecat3.hisp,shagecat4.hisp)
shagecats.asi<- c(shagecat1.asi,shagecat2.asi,shagecat3.asi,shagecat4.asi)

shagecats.races <- c(shagecats.white,shagecats.black,shagecats.hisp,shagecats.asi)

white.meanop <- sum((shagecats.races*opinion)[1:4])
white.meanop
black.meanop <- sum((shagecats.races*opinion)[5:8])
black.meanop
hisp.meanop <- sum((shagecats.races*opinion)[9:12])
hisp.meanop
asi.meanop <- sum((shagecats.races*opinion)[13:16])
asi.meanop

black.meanop.adj <- sum(shagecats.races[1:4]*opinion[5:8])
black.meanop.adj
black.meanop.adj.vec <- shagecats.races[1:4]*opinion[5:8]
hisp.meanop.adj <- sum(shagecats.races[1:4]*opinion[9:12])
hisp.meanop.adj
hisp.meanop.adj.vec <- shagecats.races[1:4]*opinion[9:12]
asi.meanop.adj <- sum(shagecats.races[1:4]*opinion[13:16])
asi.meanop.adj
asi.meanop.adj.vec <- (shagecats.races[1:4]*opinion[13:16])


black.var <- c(black.agecat1.var,black.agecat2.var,black.agecat3.var,black.agecat4.var)
hisp.var <- c(hisp.agecat1.var,hisp.agecat2.var,hisp.agecat3.var,hisp.agecat4.var)
asi.var <- c(asi.agecat1.var,asi.agecat2.var,asi.agecat3.var,asi.agecat4.var)

black.sd.adj <- sqrt(sum((shagecats.races[1:4]^2)*black.var))
black.sd.adj
hisp.sd.adj <- sqrt(sum((shagecats.races[1:4]^2)*hisp.var))
hisp.sd.adj
asi.sd.adj <- sqrt(sum((shagecats.races[1:4]^2)*asi.var))
asi.sd.adj

black.sd <- sqrt((black.meanop*(1-black.meanop))/length(incpolice[black==1]))
black.sd
white.sd <- sqrt((white.meanop*(1-white.meanop))/length(incpolice[white==1]))
white.sd
hisp.sd <- sqrt((hisp.meanop*(1-hisp.meanop))/length(incpolice[hisp==1]))
hisp.sd
asi.sd <- sqrt((asi.meanop*(1-asi.meanop))/length(incpolice[asi==1]))
asi.sd

white.sd <- sqrt((white.meanop*(1-white.meanop))/length(incpolice[white==1]))

meanop.adj <- c(black.meanop,black.meanop.adj,hisp.meanop,hisp.meanop.adj,asi.meanop,asi.meanop.adj)
races.foradj <- c("Black","Black (adjusted)","Hispanic","Hispanic (adjusted)","Asian","Asian (adjusted)")
sds.adj <- c(black.sd,black.sd.adj,hisp.sd,hisp.sd.adj,asi.sd,asi.sd.adj)

dat = data.frame(meanop.adj,races.foradj,sds.adj)
require(ggplot2)
p <- ggplot(dat,aes(races.foradj,meanop.adj,fill=races.foradj))
p +guides(fill=guide_legend(title="")) + xlab("")+ylab("")+ylim(0,.9)+
            geom_hline(yintercept = white.meanop, col = "navy", size = 1, linetype = "dashed") +
geom_bar(stat = "identity",position="dodge")+
  geom_errorbar(aes(ymin=meanop.adj-qnorm(.975)*sds.adj,ymax=meanop.adj+qnorm(.975)*sds.adj), width=.2,
                 position=position_dodge(.9))+ 
scale_fill_manual(values=c("#E69F00","#E69F00","#56B4E9","#56B4E9","#009E73","#009E73","#F0E442"))+
theme(
 panel.background = element_rect(fill = "white"), 
                 #  plot.background = element_rect(fill = "white", color = "white"), 
                  # panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                panel.grid.minor = element_blank(), 
                axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_text(size = 15), 
                 axis.text.y = element_text(size = 15))

## TABLE 2: Predicting support for additional policing

agecat <- as.factor(agecat)
racegrp <- as.factor(racegrp)
white <- as.factor(white)
black <- as.factor(black)
hisp <- as.factor(hisp)
asi <- as.factor(asi)
urban.3 <- as.factor(urban.3) 
inc <- as.factor(inc) 
ch18 <- as.factor(ch18)
own <- as.factor(own)
dem <- as.factor(dem)
ed <- as.factor(ed)
log.vicrime <- log(vicrime+1) 
log.propcrime <- log(propcrime+1) 


labs <- c("Age 30-44", "Age 45-64", "Age over 65", "Black", "Hispanic",
"Asian","Democrat","Male","Medium metro county", "Urban county",
"Income $30k-$60k","Income $60k-$100k","Income over $100k",
"HS grad","Some college", "College grad","Postgrad", "Child under 18",
"Homeowner","Violent crime rate (log)")

mainreg.col3 <- lm(incpolice~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
mainreg.col2 <- lm(incpolice~agecat+racegrp+as.factor(state),weights=comwt)
mainreg.col1 <- lm(incpolice~agecat+as.factor(state),weights=comwt)

require(stargazer)
stargazer(mainreg.col1,mainreg.col2,mainreg.col3,omit="state",no.space=T,covariate.labels=labs)

## FIGURE 3: Model-predicted responses: support for additional policing within race and age groups
require(visreg)
require(gridExtra)

p1 <- visreg(mainreg.col3,"agecat",partial=F,gg=T,xlab="Age categories",ylab="Predicted support for additional police")
p2 <- visreg(mainreg.col3,"racegrp",partial=F,gg=T,xlab="Race groups",ylab="Predicted support for additional police")
grid.arrange(p1+ylim(.2,.9),p2+ylim(.2,.9),nrow = 1)

## TABLE 3: GSS Analysis

rm(list=ls())

data <- read.csv("C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/Replication/masterdata_gss.csv")
attach(data)

summary(pred.splaw.age <- lm(splaw~birth.1898to1910+birth.1910to1922+birth.1922to1934+
birth.1934to1946+birth.1946to1958+birth.1958to1970+birth.1970to1982+ba+
dem+loginc+black+male+as.factor(region)+ch18+own+as.factor(year)+age,weights=wts))

stargazer(pred.splaw.age,no.space=T,omit=c("region","year"))

## FIGURE 4: Ratios of voter registration and political participation to age group shares of race groups

rm(list=ls()) 
data<- read.csv("C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/Replication/masterdata_cces_vv.csv")
attach(data)


shreg.age1 <- count(mactive[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(mactive,wt_var="comwt")[2,2]
shreg.age2 <- count(mactive[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(mactive,wt_var="comwt")[2,2]
shreg.age3 <- count(mactive[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(mactive,wt_var="comwt")[2,2]
shreg.age4 <- count(mactive[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(mactive,wt_var="comwt")[2,2]

shreg.age1
shreg.age2
shreg.age3
shreg.age4

shpart.age1 <- count(part[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(part,wt_var="comwt")[2,2]
shpart.age2 <- count(part[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(part,wt_var="comwt")[2,2]
shpart.age3 <- count(part[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(part,wt_var="comwt")[2,2]
shpart.age4 <- count(part[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(part,wt_var="comwt")[2,2]

shpart.age1
shpart.age2
shpart.age3
shpart.age4

shagecat1 <- weighted.mean(agecat==1,w=comwt,na.rm=T)
shagecat2 <- weighted.mean(agecat==2,w=comwt,na.rm=T)
shagecat3 <- weighted.mean(agecat==3,w=comwt,na.rm=T)
shagecat4 <- weighted.mean(agecat==4,w=comwt,na.rm=T)

shreg.agecats <- c(shreg.agecat1,shreg.agecat2,shreg.agecat3,shreg.agecat4)
shagecats <- c(shagecat1,shagecat2,shagecat3,shagecat4)

shreg.age1.white <- count(mactive[agecat==1&white==1],wt_var="comwt[agecat==1&white==1]")[2,2]/count(mactive[white==1],wt_var="comwt[white==1]")[2,2]
shreg.age1.white
shreg.age2.white <- count(mactive[agecat==2&white==1],wt_var="comwt[agecat==2&white==1]")[2,2]/count(mactive[white==1],wt_var="comwt[white==1]")[2,2]
shreg.age2.white
shreg.age3.white <- count(mactive[agecat==3&white==1],wt_var="comwt[agecat==3&white==1]")[2,2]/count(mactive[white==1],wt_var="comwt[white==1]")[2,2]
shreg.age3.white
shreg.age4.white <- count(mactive[agecat==4&white==1],wt_var="comwt[agecat==4&white==1]")[2,2]/count(mactive[white==1],wt_var="comwt[white==1]")[2,2]
shreg.age4.white

shreg.agecats.white <- c(shreg.age1.white,shreg.age2.white,shreg.age3.white,shreg.age4.white)

shreg.age1.black <- count(mactive[agecat==1&black==1],wt_var="comwt[agecat==1&black==1]")[2,2]/count(mactive[black==1],wt_var="comwt[black==1]")[2,2]
shreg.age1.black 
shreg.age2.black <- count(mactive[agecat==2&black==1],wt_var="comwt[agecat==2&black==1]")[2,2]/count(mactive[black==1],wt_var="comwt[black==1]")[2,2]
shreg.age2.black 
shreg.age3.black <- count(mactive[agecat==3&black==1],wt_var="comwt[agecat==3&black==1]")[2,2]/count(mactive[black==1],wt_var="comwt[black==1]")[2,2]
shreg.age3.black 
shreg.age4.black <- count(mactive[agecat==4&black==1],wt_var="comwt[agecat==4&black==1]")[2,2]/count(mactive[black==1],wt_var="comwt[black==1]")[2,2]
shreg.age4.black 

shreg.age1.hisp <- count(mactive[agecat==1&hisp==1],wt_var="comwt[agecat==1&hisp==1]")[2,2]/count(mactive[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shreg.age1.hisp
shreg.age2.hisp<- count(mactive[agecat==2&hisp==1],wt_var="comwt[agecat==2&hisp==1]")[2,2]/count(mactive[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shreg.age2.hisp
shreg.age3.hisp<- count(mactive[agecat==3&hisp==1],wt_var="comwt[agecat==3&hisp==1]")[2,2]/count(mactive[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shreg.age3.hisp
shreg.age4.hisp<- count(mactive[agecat==4&hisp==1],wt_var="comwt[agecat==4&hisp==1]")[2,2]/count(mactive[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shreg.age4.hisp

shreg.age1.asi <- count(mactive[agecat==1&asi==1],wt_var="comwt[agecat==1&asi==1]")[2,2]/count(mactive[asi==1],wt_var="comwt[asi==1]")[2,2]
shreg.age1.asi
shreg.age2.asi<- count(mactive[agecat==2&asi==1],wt_var="comwt[agecat==2&asi==1]")[2,2]/count(mactive[asi==1],wt_var="comwt[asi==1]")[2,2]
shreg.age2.asi
shreg.age3.asi<- count(mactive[agecat==3&asi==1],wt_var="comwt[agecat==3&asi==1]")[2,2]/count(mactive[asi==1],wt_var="comwt[asi==1]")[2,2]
shreg.age3.asi
shreg.age4.asi<- count(mactive[agecat==4&asi==1],wt_var="comwt[agecat==4&asi==1]")[2,2]/count(mactive[asi==1],wt_var="comwt[asi==1]")[2,2]
shreg.age4.asi

shreg.agecats.black <- c(shreg.age1.black,shreg.age2.black,shreg.age3.black,shreg.age4.black)
shreg.agecats.hisp <- c(shreg.age1.hisp,shreg.age2.hisp,shreg.age3.hisp,shreg.age4.hisp)
shreg.agecats.asi <- c(shreg.age1.asi,shreg.age2.asi,shreg.age3.asi,shreg.age4.asi)

require(plyr)
shagecat1.white <- count(white[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat1.white
shagecat2.white <- count(white[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat2.white
shagecat3.white <- count(white[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat3.white
shagecat4.white <- count(white[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(white,wt_var="comwt")[2,2]
shagecat4.white


shagecat1.black <- count(black[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat1.black
shagecat2.black<- count(black[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat2.black
shagecat3.black<- count(black[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat3.black
shagecat4.black<- count(black[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(black,wt_var="comwt")[2,2]
shagecat4.black

shagecat1.hisp <- count(hisp[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat1.hisp
shagecat2.hisp<- count(hisp[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat2.hisp
shagecat3.hisp<- count(hisp[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat3.hisp
shagecat4.hisp<- count(hisp[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(hisp,wt_var="comwt")[2,2]
shagecat4.hisp

shagecat1.asi <- count(asi[agecat==1],wt_var="comwt[agecat==1]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat1.asi
shagecat2.asi<- count(asi[agecat==2],wt_var="comwt[agecat==2]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat2.asi
shagecat3.asi<- count(asi[agecat==3],wt_var="comwt[agecat==3]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat3.asi
shagecat4.asi<- count(asi[agecat==4],wt_var="comwt[agecat==4]")[2,2]/count(asi,wt_var="comwt")[2,2]
shagecat4.asi



shpart.age1.white <- count(part[agecat==1&white==1],wt_var="comwt[agecat==1&white==1]")[2,2]/count(part[white==1],wt_var="comwt[white==1]")[2,2]
shpart.age1.white
shpart.age2.white <- count(part[agecat==2&white==1],wt_var="comwt[agecat==2&white==1]")[2,2]/count(part[white==1],wt_var="comwt[white==1]")[2,2]
shpart.age2.white
shpart.age3.white <- count(part[agecat==3&white==1],wt_var="comwt[agecat==3&white==1]")[2,2]/count(part[white==1],wt_var="comwt[white==1]")[2,2]
shpart.age3.white
shpart.age4.white <- count(part[agecat==4&white==1],wt_var="comwt[agecat==4&white==1]")[2,2]/count(part[white==1],wt_var="comwt[white==1]")[2,2]
shpart.age4.white

shpart.age1.black <- count(part[agecat==1&black==1],wt_var="comwt[agecat==1&black==1]")[2,2]/count(part[black==1],wt_var="comwt[black==1]")[2,2]
shpart.age1.black 
shpart.age2.black <- count(part[agecat==2&black==1],wt_var="comwt[agecat==2&black==1]")[2,2]/count(part[black==1],wt_var="comwt[black==1]")[2,2]
shpart.age2.black 
shpart.age3.black <- count(part[agecat==3&black==1],wt_var="comwt[agecat==3&black==1]")[2,2]/count(part[black==1],wt_var="comwt[black==1]")[2,2]
shpart.age3.black 
shpart.age4.black <- count(part[agecat==4&black==1],wt_var="comwt[agecat==4&black==1]")[2,2]/count(part[black==1],wt_var="comwt[black==1]")[2,2]
shpart.age4.black 

shpart.age1.hisp <- count(part[agecat==1&hisp==1],wt_var="comwt[agecat==1&hisp==1]")[2,2]/count(part[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shpart.age1.hisp
shpart.age2.hisp<- count(part[agecat==2&hisp==1],wt_var="comwt[agecat==2&hisp==1]")[2,2]/count(part[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shpart.age2.hisp
shpart.age3.hisp<- count(part[agecat==3&hisp==1],wt_var="comwt[agecat==3&hisp==1]")[2,2]/count(part[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shpart.age3.hisp
shpart.age4.hisp<- count(part[agecat==4&hisp==1],wt_var="comwt[agecat==4&hisp==1]")[2,2]/count(part[hisp==1],wt_var="comwt[hisp==1]")[2,2]
shpart.age4.hisp

shpart.age1.asi <- count(part[agecat==1&asi==1],wt_var="comwt[agecat==1&asi==1]")[2,2]/count(part[asi==1],wt_var="comwt[asi==1]")[2,2]
shpart.age1.asi
shpart.age2.asi<- count(part[agecat==2&asi==1],wt_var="comwt[agecat==2&asi==1]")[2,2]/count(part[asi==1],wt_var="comwt[asi==1]")[2,2]
shpart.age2.asi
shpart.age3.asi<- count(part[agecat==3&asi==1],wt_var="comwt[agecat==3&asi==1]")[2,2]/count(part[asi==1],wt_var="comwt[asi==1]")[2,2]
shpart.age3.asi
shpart.age4.asi<- count(part[agecat==4&asi==1],wt_var="comwt[agecat==4&asi==1]")[2,2]/count(part[asi==1],wt_var="comwt[asi==1]")[2,2]
shpart.age4.asi

shagecats.white <- c(shagecat1.white,shagecat2.white,shagecat3.white,shagecat4.white)
shagecats.black <- c(shagecat1.black,shagecat2.black,shagecat3.black,shagecat4.black)
shagecats.hisp<- c(shagecat1.hisp,shagecat2.hisp,shagecat3.hisp,shagecat4.hisp)
shagecats.asi<- c(shagecat1.asi,shagecat2.asi,shagecat3.asi,shagecat4.asi)

shagecats.races <- c(shagecats.white,shagecats.black,shagecats.hisp,shagecats.asi)

shpart.agecats.white <- c(shpart.age1.white,shpart.age2.white,shpart.age3.white,shpart.age4.white)
shpart.agecats.black <- c(shpart.age1.black,shpart.age2.black,shpart.age3.black,shpart.age4.black)
shpart.agecats.hisp <- c(shpart.age1.hisp,shpart.age2.hisp,shpart.age3.hisp,shpart.age4.hisp)
shpart.agecats.asi <- c(shpart.age1.asi,shpart.age2.asi,shpart.age3.asi,shpart.age4.asi)

shpart.agecats.races <- c(shpart.agecats.white,shpart.agecats.black,shpart.agecats.hisp,shpart.agecats.asi)
shreg.agecats.races <- c(shreg.agecats.white,shreg.agecats.black,shreg.agecats.hisp,shreg.agecats.asi)

Races <-c(rep("White", 4), rep("Black", 4), rep("Hispanic",4), rep("Asian",4))
Races
agecats <- c("<30","30-44","45-64","65+","<30","30-44","45-64","65+","<30","30-44","45-64","65+","<30","30-44","45-64","65+")
agecats <- as.factor(agecats)

partcomp <- data.frame(races,shagecats.races,shreg.agecats.races,shpart.agecats.races)


partgraph.write <- write.csv(partcomp,"C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/partcomp.csv")

regratio.white <- shreg.agecats.white/shagecats.white
regratio.black<- shreg.agecats.black/shagecats.black
regratio.hisp<- shreg.agecats.hisp/shagecats.hisp
regratio.asi<- shreg.agecats.asi/shagecats.asi

partratio.white <- shpart.agecats.white/shagecats.white
partratio.black <- shpart.agecats.black/shagecats.black
partratio.hisp <- shpart.agecats.hisp/shagecats.hisp
partratio.asi <- shpart.agecats.asi/shagecats.asi

regratios <- c(regratio.white,regratio.black,regratio.hisp,regratio.asi)
regratios
partratios <- c(partratio.white,partratio.black,partratio.hisp,partratio.asi)
partratios 

dat <- data.frame(regratios,Races,agecats)
regratios.write <- write.csv(dat,"C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/regratios.csv")

colnames(dat)[colnames(dat)=="races"] <- "Race groups"

dat2 <- data.frame(partratios,Races,agecats)


p <-ggplot(dat, aes(agecats,regratios))
p +guides(fill=guide_legend(title=""))+xlab("Age groups") + ylab("Registration ratio")+ylim(0,2)+geom_hline(yintercept = 1, col = "darkblue", size = 1, linetype = "solid")+
geom_bar(stat = "identity", aes(fill = Races),position="dodge")+theme(
 panel.background = element_rect(fill = "white"), 
                 #  plot.background = element_rect(fill = "white", color = "white"), 
                  # panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                panel.grid.minor = element_blank(), 
                axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_text(size = 15), 
                 axis.text.y = element_text(size = 15))

q <- ggplot(dat2,aes(agecats,partratios))
q +guides(fill=guide_legend(title=""))+xlab("Age groups") + ylab("Other participation ratio")+ylim(0,2)+geom_hline(yintercept = 1, col = "darkblue", size = 1, linetype = "solid")+
geom_bar(stat = "identity", aes(fill = Races),position="dodge")+theme(
 panel.background = element_rect(fill = "white"), 
                 #  plot.background = element_rect(fill = "white", color = "white"), 
                  # panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                panel.grid.minor = element_blank(), 
                axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_text(size = 15), 
                 axis.text.y = element_text(size = 15))


grid.arrange(p +ggtitle("Registration ratios")+guides(fill=guide_legend(title=""))+xlab("") + ylab("Registration ratio")+ylim(0,2)+
geom_bar(stat = "identity", aes(fill = Races),position="dodge")+geom_hline(yintercept = 1, col = "darkblue", size = 1, linetype = "solid")+theme(
 panel.background = element_rect(fill = "white"), 
                 #  plot.background = element_rect(fill = "white", color = "white"), 
                  # panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                panel.grid.minor = element_blank(), 
                axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_blank(), 
                 axis.text.y = element_text(size = 15)),
q+ ggtitle("Other participation ratios")+guides(fill=guide_legend(title=""))+xlab("Age groups") + ylab("Other participation ratio")+ylim(0,2)+
geom_bar(stat = "identity", aes(fill = Races),position="dodge")+geom_hline(yintercept = 1, col = "darkblue", size = 1, linetype = "solid")+theme(
 panel.background = element_rect(fill = "white"), 
                 #  plot.background = element_rect(fill = "white", color = "white"), 
                  # panel.grid.major = element_blank(), 
                  panel.border = element_blank(), 
                panel.grid.minor = element_blank(), 
                axis.line.x = element_line(color = "gray40"),
                  axis.text.x = element_text(size = 15), 
                  axis.title.y = element_blank(), 
                 axis.text.y = element_text(size = 15)))

## LOOK HERE - CHANGE THIS IN THE DISSERTATION AND PAPER
## TABLE 4: Predicting registered voter status and non-voting political participation

agecat <- as.factor(agecat)
racegrp <- as.factor(racegrp)
white <- as.factor(white)
black <- as.factor(black)
hisp <- as.factor(hisp)
asi <- as.factor(asi)
urban.3 <- as.factor(urban.3) 
inc <- as.factor(inc) 
child18 <- as.factor(child18)
own <- as.factor(own)
dem <- as.factor(dem)
ed <- as.factor(ed)
log.vicrime <- log(vicrime+1) 
log.propcrime <- log(propcrime+1) 

pred.mactive <- lm(mactive~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+child18+own+log.vicrime+as.factor(state),weights=comwt)
summary(pred.mactive)

pred.part <- lm(part~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+child18+own+log.vicrime+as.factor(state),weights=comwt)
summary(pred.part)

stargazer(pred.mactive,pred.part,omit="state",no.space=T)

########## Cities analysis

## TABLE 5: Summary stats for 1,539 cities
rm(list=ls())

data <- read.csv("C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/Replication/masterdata_cities.csv")
attach(data)

stargazer(data)

## TABLE 6: Predicting police force size

log.shover65 <- log(shover65)
log.pcinc <- log(pcinc)
log.totalpop <- log(totalpop)
log.vicrime <- log(vicrime+1)

summary(fit.logopp <- lm(log.opp~log.shover65+shwhite+log.pcinc+log.vicrime+log.totalpop+as.factor(state)))

stargazer(fit.logopp,omit="state")

## FIGURE 5: Predictive exercise: full-time sworn ocers per 1,000 persons

p1 <- visreg(fit.logopp,partial=F,ylim=c(0,10),xvar="log.shover65",gg=F,trans=exp,xtrans=exp,
xlab="Share over 65", ylab="Officers per person")
plot(p1,ylim=c(0,10),xlab="Share over 65",ylab="FT sworn officers per 1000 persons")

## REPLICATING APPENDIX TABLES

## A1: Replication of Table 1 w/ continuous age

rm(list=ls())

data <- read.csv("C:/Users/Becca/Dropbox/CCES 2016 - age cohort differences/Data/Replication/masterdata_cces.csv")
attach(data)

agecat <- as.factor(agecat)
racegrp <- as.factor(racegrp)
white <- as.factor(white)
black <- as.factor(black)
hisp <- as.factor(hisp)
asi <- as.factor(asi)
urban.3 <- as.factor(urban.3) 
inc <- as.factor(inc) 
ch18 <- as.factor(ch18)
own <- as.factor(own)
dem <- as.factor(dem)
ed <- as.factor(ed)
log.vicrime <- log(vicrime+1) 
log.propcrime <- log(propcrime+1) 

labs.continuous <- c("Age (years)", "Black", "Hispanic",
"Asian","Democrat","Male","Medium metro county", "Urban county",
"Income $30k-$60k","Income $60k-$100k","Income over $100k",
"HS grad","Some college", "College grad","Postgrad", "Child under 18",
"Homeowner","Violent crime rate (log)")
mainreg.age.continuous <- lm(incpolice~age+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
stargazer(mainreg.age.continuous,omit="state",no.space=T,covariate.labels=labs.continuous)

## A2: Replication of Table 1 w/ robust standard errors

library(lmtest)
library(sandwich)

mainreg.col3 <- lm(incpolice~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
mainreg.col2 <- lm(incpolice~agecat+racegrp+as.factor(state),weights=comwt)
mainreg.col1 <- lm(incpolice~agecat+as.factor(state),weights=comwt)

robust.ses.1 <- coeftest(mainreg.col1, vcov = vcovHC(mainreg.col1, type="HC1"))
robust.ses.1
robust.ses.1[,2:3]

robust.ses.2 <- coeftest(mainreg.col2, vcov = vcovHC(mainreg.col2, type="HC1"))
robust.ses.2
robust.ses.2[,2:3]

robust.ses.3 <- coeftest(mainreg.col3, vcov = vcovHC(mainreg.col3, type="HC1"))
robust.ses.3
robust.ses.3[,2:3]

## A3: Predict police feeling, support for increased education spending, traffic ticket receipt, and recent victimization

pred.policefeel <- lm(policefeel~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
summary(pred.policefeel)

pred.edspend.bin<- lm(edspend.bin~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
summary(pred.edspend.bin)

pred.ticket<- lm(ticket~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
summary(pred.ticket)

pred.victim <- lm(victim~agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
summary(pred.victim)

labs <- c("Age 30-44", "Age 45-64", "Age over 65", "Black", "Hispanic",
"Asian","Democrat","Male","Medium metro county", "Urban county",
"Income $30k-$60k","Income $60k-$100k","Income over $100k",
"HS grad","Some college", "College grad","Postgrad", "Child under 18",
"Homeowner","Violent crime rate (log)")


stargazer(pred.policefeel,pred.edspend.bin,pred.ticket,pred.victim,omit="state",no.space=T,covariate.labels=labs)

## A4: Predicting support for additional policing with these four variables

labs.meds <- c("Police make R feel safe", "Support for increased education spending", "Traffic ticket receipt (last 2 yrs)", "Crime victim (last 2 yrs)", labs)
mainreg.policefeel <- lm(incpolice~policefeel+agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
mainreg.edspend <- lm(incpolice~edspend.bin+agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
mainreg.ticket <- lm(incpolice~ticket+agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)
mainreg.vic <- lm(incpolice~victim+agecat+racegrp+dem+male+urban.3+(inc)+(ed)+ch18+own+log.vicrime+as.factor(state),weights=comwt)

stargazer(mainreg.policefeel,mainreg.edspend,mainreg.ticket,mainreg.vic,omit="state",no.space=T,covariate.labels=labs.meds)


